#include <stdio.h>
#include <string.h>
#include <math.h>



    /* Funcao que executa o algoritmo de Gauss_Jacobi */
__global__ void Gauss_Jacobi_GPU(int* iter,int iter_max, float* A, float* x, float* b, int order, int error) {
    int index = blockDim.x * blockIdx.x + threadIdx.x; 
    
    if (index < order) {
        __shared__ float numerator;
        __shared__ float denum;
        numerator = 1;
        denum = 1;
        float x1;
        iter[index] = 0;
        while ((numerator / denum) > error && iter[index] < iter_max) {
            numerator = 0;
            denum = 0;
            x1 = 0;
            iter[index]++;
            for (int k = 0; k < order; ++k) {
                if (index != k)
                    x1 += A[(index * order) + k] * x[k];
            }
            
            x1 = (b[index] - x1)/ A[(index * order) + index];

            if (numerator < abs(x1 - x[index]))
                numerator = abs(x1 - x[index]);
            if (denum < abs(x1))
                denum = abs(x1);
                
            x[index] = x1;
            __syncthreads(); /* Para evitar erro por divisão em alguns casos */
        }
    }
}

    /* Funcao que efetua o teste da convergencia na matriz A */
__global__ void ConvergenceLine_Test_GPU(int *Convergence,float* A, float* x, float* b, int order){

    int index = blockDim.x * blockIdx.x + threadIdx.x; 

    if (index < order) {
        __shared__ float maxl;
        __shared__ float suml;
        maxl = 0;
        int j;

        suml = 0;
        for(j = 0; j < order; j++)
            if(index != j)
                suml += A[(index * order) + j];
        
        suml = suml / abs(A[index*order + index]);
        if(maxl < suml)
            maxl = suml;

          if(maxl > 1)
              *Convergence = 0;
          else
              *Convergence = 1;
    }
}




int main(int argc, char** argv) {

FILE *ifp; /* Variavel que contem o ponteiro para arquivo de entrada */
    
    /* Verifica se foi possivel abrir o arquivo */
    if((ifp = fopen(argv[1], "r")) == NULL)
    {
        printf("Could not open the file %s\n", argv[1]);
    }

    int order = 0; /* Variavel que indica a ordem da matriz */
    int row_test = 0; /* Variavel que indica a linha a ser avaliada */
    int i; /* Variavel contadora */
    int j; /* Variavel contadora */
    int ite_max = 0; /* Variavel que indica o numero maximo de iteracoes */

    float error = 0.0f; /* Variavel que indica erro para a condicao de parada */

     /* Recebe valores de definicao da matriz */
    fscanf(ifp, "%d\n%d\n%f\n%d", &order, &row_test, &error, &ite_max);
    
    float* Matrix_A = (float*) malloc(order*order*sizeof(float)); /* Declaracao da matriz A usada na CPU */
    float* Vector_x = (float*) malloc(order*sizeof(float)); /* Declaracao do vetor x usado na CPU */
    float* Vector_b = (float*) malloc(order*sizeof(float)); /* Declaracao do vetor b usado na CPU */
    float* Result_CPU = (float*) malloc(order*sizeof(float)); /* Declaracao da vetor resultado na CPU */
    int*   Iter_CPU = (int*) malloc(order*sizeof(int)); /* Declaracao usado para saber as iteracoes */

    /* Declaracao da Matriz e vetores utilizados na placa */
    float* Matrix_A_Device;
    cudaMalloc((void**) &Matrix_A_Device, order*order*sizeof(float));
    float* Vector_x_device;
    cudaMalloc((void**) &Vector_x_device,order*sizeof(float));
    float* Vector_b_device;
    cudaMalloc((void**) &Vector_b_device,order*sizeof(float));
    int* Iterations_device;
    cudaMalloc((void**) &Iterations_device,order*sizeof(int));
   
   
    
    /* Recebe os valores do arquivo para a matriz A */
    for(i = 0; i < order; i++)
        for(j = 0; j < order; j++)
            fscanf(ifp, "%f", &Matrix_A[(i * order) + j]);

    /* Recebe os valores do arquivo para a matriz B */
    for(i = 0; i < order; i++)
        fscanf(ifp, "%f", &Vector_b[i]);

    /*Inicializa vetor x*/    
    for(i = 0; i < order; i++)
        Vector_x[i] = Vector_b[i]/Matrix_A[i*order + i];
        
    
    /* Copia para a placa o conteudo da Matriz A e vetores x e b para a placa*/
    cudaMemcpy(Matrix_A_Device, Matrix_A, order*order*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(Vector_x_device, Vector_x, order*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(Vector_b_device, Vector_b, order*sizeof(float), cudaMemcpyHostToDevice);

    int blocos=32; /* Define o numero de blocos utilizado*/
    
    /* Declaracao das variaveis utilizada na CPU e na GPU sobre a convergencia */
    int *PtrConvergence;
    int Convergence;
    cudaMalloc((void**) &PtrConvergence,sizeof(int));
    ConvergenceLine_Test_GPU<<< blocos, order/32 +1 >>>(PtrConvergence,Matrix_A_Device, Vector_x_device, Vector_b_device,order);
    cudaMemcpy(&Convergence,PtrConvergence, sizeof(int),cudaMemcpyDeviceToHost);

    /* Verifica o criterio de convergencia*/
    if(Convergence == 0){
      fprintf(stderr,"Convergence criterial not satisfied");
      exit(-1);
    }
    else printf("Convergence criterion satisfied.\n");
       
   
    /* Efetua o algoritmo de Gauss_Jacobi */
    Gauss_Jacobi_GPU<<< blocos, order/32 +1 >>>(Iterations_device,ite_max, Matrix_A_Device, Vector_x_device, Vector_b_device,order,error);
    
    
    /* Copia as os resultados e iteracoes para a CPU */  
    cudaMemcpy(Result_CPU, Vector_x_device, order*sizeof(float),cudaMemcpyDeviceToHost);
    cudaMemcpy(Iter_CPU, Iterations_device, order*sizeof(float),cudaMemcpyDeviceToHost);

    
    int h; /* Variavel auxiliar*/
        
    /* Verifica a maior iteraçao */
    h = -1;
    for(i = 0; i < order;i++)
      if(h < Iter_CPU[i])
        h = Iter_CPU[i];
      
    
    /* Calcula o resultado */  
    float res=0;  
    for(i = 0; i < order; i++)
        res += Matrix_A[(row_test * order) + i] * Result_CPU[i];
        
    
    /* Imprime na tela os valores do resultado e iteraçao */    
    printf("-----------------------------------------\n");
    printf("Ite: %d\n",h);
    printf("RowTest: %d => [%f] =? [%f]\n",
            row_test, res, Vector_b[row_test]);
    printf("-----------------------------------------\n");

    
    /* Libera a memoria */
    free(Matrix_A);
    free(Vector_b);
    free(Result_CPU);
    free(Vector_x);
    cudaFree(PtrConvergence);
    cudaFree(Matrix_A_Device);
    cudaFree(Vector_b_device);
    cudaFree(Vector_x_device);
    cudaFree(Iterations_device);

        
    return 0;
}
